!> \file GFS_surface_generic_post.F90
!!  Contains code related to all GFS surface schemes to be run afterward.

      module GFS_surface_generic_post

      use machine, only: kind_phys

      implicit none

      private

      public GFS_surface_generic_post_init, GFS_surface_generic_post_run

      real(kind=kind_phys), parameter :: zero = 0.0_kind_phys, one = 1.0_kind_phys

      contains

!>\defgroup gfs_sfc_gen_post_mode GFS surface_generic_post Module
!! This module contains code related to all GFS surface schemes to be run afterward.
!> @{
!> \section arg_table_GFS_surface_generic_post_init Argument Table
!! \htmlinclude GFS_surface_generic_post_init.html
!!
      subroutine GFS_surface_generic_post_init (vtype, stype, slope, vtype_save, stype_save, slope_save, errmsg, errflg)

        integer, dimension(:), intent(in)  :: vtype_save, stype_save, slope_save
        integer, dimension(:), intent(out) :: vtype, stype, slope

        ! CCPP error handling
        character(len=*), intent(out) :: errmsg
        integer,          intent(out) :: errflg

        ! Initialize CCPP error handling variables
        errmsg = ''
        errflg = 0

        ! Restore vegetation, soil and slope type
        vtype(:) = vtype_save(:)
        stype(:) = stype_save(:)
        slope(:) = slope_save(:)

      end subroutine GFS_surface_generic_post_init

!> \section arg_table_GFS_surface_generic_post_run Argument Table
!! \htmlinclude GFS_surface_generic_post_run.html
!!
      subroutine GFS_surface_generic_post_run (im, cplflx, cplaqm, cplchm, cplwav, cpllnd, lssav, dry, icy, wet,                    &
        lsm, lsm_noahmp, dtf, ep1d, gflx, tgrs_1, qgrs_1, ugrs_1, vgrs_1,                                                           &
        adjsfcdlw, adjsfcdsw, adjnirbmd, adjnirdfd, adjvisbmd, adjvisdfd, adjsfculw, adjsfculw_wat, adjnirbmu, adjnirdfu,           &
        adjvisbmu, adjvisdfu, t2m, q2m, u10m, v10m, tsfc, tsfc_wat, pgr, xcosz, evbs, evcw, trans, sbsno, snowc, snohf, pah, pahi,  &
        epi, gfluxi, t1, q1, u1, v1, dlwsfci_cpl, dswsfci_cpl, dlwsfc_cpl, dswsfc_cpl, dnirbmi_cpl, dnirdfi_cpl, dvisbmi_cpl,       &
        dvisdfi_cpl, dnirbm_cpl, dnirdf_cpl, dvisbm_cpl, dvisdf_cpl, nlwsfci_cpl, nlwsfc_cpl, t2mi_cpl, q2mi_cpl, u10mi_cpl,        &
        v10mi_cpl, tsfci_cpl, psurfi_cpl, nnirbmi_cpl, nnirdfi_cpl, nvisbmi_cpl, nvisdfi_cpl, nswsfci_cpl, nswsfc_cpl, nnirbm_cpl,  &
        nnirdf_cpl, nvisbm_cpl, nvisdf_cpl, gflux, evbsa, evcwa, transa, sbsnoa, snowca, snohfa, paha, ep, ecan, etran, edir, waxy, &
        runoff, srunoff, runof, drain, tecan, tetran, tedir, twa, lheatstrg, h0facu, h0facs, zvfun, hflx, evap, hflxq, hffac,       &
        isot, ivegsrc, islmsk, vtype, stype, slope, vtype_save, stype_save, slope_save, errmsg, errflg)

        implicit none

        integer,                                intent(in) :: im
        logical,                                intent(in) :: cplflx, cplaqm, cplchm, cplwav, cpllnd, lssav
        logical, dimension(:),                  intent(in) :: dry, icy, wet
        integer,                                intent(in) :: lsm, lsm_noahmp
        real(kind=kind_phys),                   intent(in) :: dtf

        real(kind=kind_phys), dimension(:),  intent(in)  :: ep1d, gflx, tgrs_1, qgrs_1, ugrs_1, vgrs_1, adjsfcdlw, adjsfcdsw,  &
          adjnirbmd, adjnirdfd, adjvisbmd, adjvisdfd, adjsfculw, adjsfculw_wat, adjnirbmu, adjnirdfu, adjvisbmu, adjvisdfu,    &
          t2m, q2m, u10m, v10m, tsfc, tsfc_wat, pgr, xcosz, evbs, evcw, trans, sbsno, snowc, snohf, pah, ecan, etran, edir,    &
          waxy

        real(kind=kind_phys), dimension(:),  intent(inout) :: epi, gfluxi, t1, q1, u1, v1, dlwsfci_cpl, dswsfci_cpl, dlwsfc_cpl, &
          dswsfc_cpl, dnirbmi_cpl, dnirdfi_cpl, dvisbmi_cpl, dvisdfi_cpl, dnirbm_cpl, dnirdf_cpl, dvisbm_cpl, dvisdf_cpl,        &
          nlwsfci_cpl, nlwsfc_cpl, t2mi_cpl, q2mi_cpl, u10mi_cpl, v10mi_cpl, tsfci_cpl, psurfi_cpl, nnirbmi_cpl, nnirdfi_cpl,    &
          nvisbmi_cpl, nvisdfi_cpl, nswsfci_cpl, nswsfc_cpl, nnirbm_cpl, nnirdf_cpl, nvisbm_cpl, nvisdf_cpl, gflux, evbsa,       &
          evcwa, transa, sbsnoa, snowca, snohfa, ep, paha, tecan, tetran, tedir, twa, pahi

        real(kind=kind_phys), dimension(:), intent(inout) :: runoff, srunoff
        real(kind=kind_phys), dimension(:), intent(in)    :: drain, runof

        ! For canopy heat storage
        logical, intent(in) :: lheatstrg
        real(kind=kind_phys), intent(in) :: h0facu, h0facs
        real(kind=kind_phys), dimension(:), intent(in)  :: zvfun
        real(kind=kind_phys), dimension(:), intent(in)  :: hflx,  evap
        real(kind=kind_phys), dimension(:), intent(out) :: hflxq
        real(kind=kind_phys), dimension(:), intent(out) :: hffac

        integer, intent(in) :: isot, ivegsrc, islmsk(:), vtype_save(:), stype_save(:), slope_save(:)
        integer, intent(out) :: vtype(:), stype(:), slope(:)

        ! CCPP error handling variables
        character(len=*), intent(out) :: errmsg
        integer,          intent(out) :: errflg

        ! Local variables
        real(kind=kind_phys), parameter :: albdf = 0.06_kind_phys

        integer :: i
        real(kind=kind_phys) :: xcosz_loc, ocalnirdf_cpl, ocalnirbm_cpl, ocalvisdf_cpl, ocalvisbm_cpl

        ! Initialize CCPP error handling variables
        errmsg = ''
        errflg = 0

        do i=1,im
          epi(i)    = ep1d(i)
          gfluxi(i) = gflx(i)
          if (lsm == lsm_noahmp) then
            pahi(i)   = pah(i)
          endif
          t1(i)     = tgrs_1(i)
          q1(i)     = qgrs_1(i)
          u1(i)     = ugrs_1(i)
          v1(i)     = vgrs_1(i)
        enddo

        if (cplflx .or. cplchm .or. cplwav) then
          do i=1,im
            u10mi_cpl(i) = u10m(i)
            v10mi_cpl(i) = v10m(i)
          enddo
        endif

        if (cplflx .or. cplchm .or. cpllnd) then
          do i=1,im
            tsfci_cpl(i) = tsfc(i)
          enddo
        endif

        if (cplflx .or. cpllnd) then
          do i=1,im
            dlwsfc_cpl  (i) = dlwsfc_cpl(i) + adjsfcdlw(i)*dtf
            dswsfc_cpl  (i) = dswsfc_cpl(i) + adjsfcdsw(i)*dtf
            psurfi_cpl  (i) = pgr(i)
          enddo
        endif

        if (cplflx) then
          do i=1,im
            dlwsfci_cpl (i) = adjsfcdlw(i)
            dswsfci_cpl (i) = adjsfcdsw(i)
            dnirbmi_cpl (i) = adjnirbmd(i)
            dnirdfi_cpl (i) = adjnirdfd(i)
            dvisbmi_cpl (i) = adjvisbmd(i)
            dvisdfi_cpl (i) = adjvisdfd(i)
            dnirbm_cpl  (i) = dnirbm_cpl(i) + adjnirbmd(i)*dtf
            dnirdf_cpl  (i) = dnirdf_cpl(i) + adjnirdfd(i)*dtf
            dvisbm_cpl  (i) = dvisbm_cpl(i) + adjvisbmd(i)*dtf
            dvisdf_cpl  (i) = dvisdf_cpl(i) + adjvisdfd(i)*dtf
            nlwsfci_cpl (i) = adjsfcdlw(i)  - adjsfculw(i)
            if (wet(i)) then
              nlwsfci_cpl(i) = adjsfcdlw(i) - adjsfculw_wat(i)
            endif
            nlwsfc_cpl  (i) = nlwsfc_cpl(i) + nlwsfci_cpl(i)*dtf
            t2mi_cpl    (i) = t2m(i)
            q2mi_cpl    (i) = q2m(i)
          enddo
        endif

!  ---  estimate mean albedo for ocean point without ice cover and apply
!       them to net SW heat fluxes

        if (cplflx .or. cpllnd) then
          do i=1,im
!           if (Sfcprop%landfrac(i) < one) then ! Not 100% land
            if (wet(i)) then                    ! some open water
!  ---  compute open water albedo
              xcosz_loc = max( zero, min( one, xcosz(i) ))
              ocalnirdf_cpl = 0.06_kind_phys
              ocalnirbm_cpl = max(albdf, 0.026_kind_phys/(xcosz_loc**1.7_kind_phys+0.065_kind_phys)     &
       &                       + 0.15_kind_phys * (xcosz_loc-0.1_kind_phys) * (xcosz_loc-0.5_kind_phys) &
       &                       * (xcosz_loc-one))
              ocalvisdf_cpl = 0.06_kind_phys
              ocalvisbm_cpl = ocalnirbm_cpl

              nnirbmi_cpl(i) = adjnirbmd(i) * (one-ocalnirbm_cpl)
              nnirdfi_cpl(i) = adjnirdfd(i) * (one-ocalnirdf_cpl)
              nvisbmi_cpl(i) = adjvisbmd(i) * (one-ocalvisbm_cpl)
              nvisdfi_cpl(i) = adjvisdfd(i) * (one-ocalvisdf_cpl)
            else
              nnirbmi_cpl(i) = adjnirbmd(i) - adjnirbmu(i)
              nnirdfi_cpl(i) = adjnirdfd(i) - adjnirdfu(i)
              nvisbmi_cpl(i) = adjvisbmd(i) - adjvisbmu(i)
              nvisdfi_cpl(i) = adjvisdfd(i) - adjvisdfu(i)
            endif
            nswsfci_cpl(i) = nnirbmi_cpl(i) + nnirdfi_cpl(i)   &
                           + nvisbmi_cpl(i) + nvisdfi_cpl(i)
            nswsfc_cpl(i)  = nswsfc_cpl(i)  + nswsfci_cpl(i)*dtf
            nnirbm_cpl(i)  = nnirbm_cpl(i)  + nnirbmi_cpl(i)*dtf
            nnirdf_cpl(i)  = nnirdf_cpl(i)  + nnirdfi_cpl(i)*dtf
            nvisbm_cpl(i)  = nvisbm_cpl(i)  + nvisbmi_cpl(i)*dtf
            nvisdf_cpl(i)  = nvisdf_cpl(i)  + nvisdfi_cpl(i)*dtf
          enddo
        endif

        if (cplaqm .and. .not.cplflx) then
          do i=1,im
            t2mi_cpl    (i) = t2m(i)
            q2mi_cpl    (i) = q2m(i)
            psurfi_cpl  (i) = pgr(i)
            if (wet(i)) then                    ! some open water
!  ---  compute open water albedo
              xcosz_loc = max( zero, min( one, xcosz(i) ))
              ocalnirdf_cpl = 0.06_kind_phys
              ocalnirbm_cpl = max(albdf, 0.026_kind_phys/(xcosz_loc**1.7_kind_phys+0.065_kind_phys)     &
       &                       + 0.15_kind_phys * (xcosz_loc-0.1_kind_phys) * (xcosz_loc-0.5_kind_phys) &
       &                       * (xcosz_loc-one))
              ocalvisdf_cpl = 0.06_kind_phys
              ocalvisbm_cpl = ocalnirbm_cpl

              nswsfci_cpl(i) = adjnirbmd(i) * (one-ocalnirbm_cpl) + &
                               adjnirdfd(i) * (one-ocalnirdf_cpl) + &
                               adjvisbmd(i) * (one-ocalvisbm_cpl) + &
                               adjvisdfd(i) * (one-ocalvisdf_cpl)
            else
              nswsfci_cpl(i) = adjnirbmd(i) - adjnirbmu(i) + &
                               adjnirdfd(i) - adjnirdfu(i) + &
                               adjvisbmd(i) - adjvisbmu(i) + &
                               adjvisdfd(i) - adjvisdfu(i)
            endif
          enddo
        endif

        if (lssav) then
          do i=1,im
            gflux(i)   = gflux(i)  + gflx(i)  * dtf
            evbsa(i)   = evbsa(i)  + evbs(i)  * dtf
            evcwa(i)   = evcwa(i)  + evcw(i)  * dtf
            transa(i)  = transa(i) + trans(i) * dtf
            sbsnoa(i)  = sbsnoa(i) + sbsno(i) * dtf
            snowca(i)  = snowca(i) + snowc(i) * dtf
            snohfa(i)  = snohfa(i) + snohf(i) * dtf
            ep(i)      = ep(i)     + ep1d(i)  * dtf

!  --- ...  total runoff is composed of drainage into water table and
!           runoff at the surface and is accumulated in unit of meters
            runoff(i)  = runoff(i)  + (drain(i)+runof(i)) * dtf
            srunoff(i) = srunoff(i) + runof(i) * dtf
            tecan(i)   = tecan(i)   + ecan(i) * dtf
            tetran(i)  = tetran(i)  + etran(i) * dtf
            tedir(i)   = tedir(i)   + edir(i) * dtf
            if (lsm == lsm_noahmp) then
             paha(i)    = paha(i)    + pah(i)   * dtf
             twa(i)     = waxy(i) 
            endif
          enddo
        endif

!
!  in order to achieve heat storage within canopy layer, in the canopy
!    heat torage parameterization the kinematic sensible heat flux
!    (hflx) as surface boundary forcing to the pbl scheme is
!    reduced in a factor of hffac given as a function of surface roughness &
!    green vegetation fraction (zvfun) 
!
        do i=1,im
          hflxq(i) = hflx(i)
          hffac(i) = 1.0
        enddo
        if (lheatstrg) then
          do i=1,im
            if (dry(i)) then
              if(hflx(i) > 0.) then
                hffac(i) = h0facu * zvfun(i)
              else
                hffac(i) = h0facs * zvfun(i)
              endif
              hffac(i) = 1. + hffac(i)
              hflxq(i) = hflx(i) / hffac(i)
            endif
          enddo
        endif

        ! Restore vegetation, soil and slope type
        vtype(:) = vtype_save(:)
        stype(:) = stype_save(:)
        slope(:) = slope_save(:)

      end subroutine GFS_surface_generic_post_run
!> @}
      end module GFS_surface_generic_post
